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Abstract 

A conservative discretization of incompressible Navier-Stokes equations is de¬ 
veloped based on discrete exterior calculus (DEC). A distinguishing feature of 
our method is the use of an algebraic discretization of the interior product op¬ 
erator and a combinatorial discretization of the wedge product. The governing 
equations are first rewritten using the exterior calculus notation, replacing vector 
calculus differential operators by the exterior derivative, Hodge star and wedge 
product operators. The discretization is then carried out by substituting with the 
corresponding discrete operators based on the DEC framework. Numerical exper¬ 
iments for flows over surfaces reveal a second order accuracy for the developed 
scheme when using structured-triangular meshes, and first order accuracy for oth¬ 
erwise unstructured meshes. By construction, the method is conservative in that 
both mass and vorticity are conserved up to machine precision. The relative error 
in kinetic energy for inviscid flow test cases converges in a second order fashion 
with both the mesh size and the time step. 
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1. Introduction 


When solving a partial differential equation numerically, various measures 
(e.g. convergence, stability and consistency) are usually investigated to verify 
the implemented discretization. Such measures, although reflecting the mathe¬ 
matical fidelity of the discretization, may not give insight into the physical fidelity 
of the discretization. By physical fidelity we mean how well does the discrete sys¬ 
tem of equations conserve secondary quantities, such as kinetic energy, implied 
in the continuous equation but not explicitly constructed or built into the numeri¬ 
cal scheme. The development of such physically conservative discretizations, for 
Navier-Stokes (NS) equations for example, is favorable for many physical appli¬ 
cations (e.g. turbulent flows and shallow-water simulations) to avoid undesirable 
numerical artifacts. Among other discretization approaches, some staggered mesh 
schemes are known for their conservation of both primary (i.e. mass and momen¬ 
tum) and secondary (e.g. kinetic energy and vorticity) physical quantities (B. 

Staggered mesh methods were first developed by Harlow and Welch |2l for 
structured Cartesian meshes by placing the velocity and pressure degrees of free¬ 
dom at different positions on the mesh. Later on, the approach was extended to 
unstructured meshes by Nicolaides 0 and Hall et al. 01, which is now known as 
the covolume (or dual-variable) method. The covolume method was originally in¬ 
troduced as a low order method that does not experience spurious modes that were 
common in low order discretizations of viscous flows. The derivation of such dis¬ 
cretization commences by taking the dot product of the momentum equation with 
the unit normal vector to each face of the triangular/tetrahedral elements. This 
reduces the velocity vector field to a scalar flux (equal to the mass flux across the 
face for an incompressible flow with constant density) defined on each face. In 
this approach, the pressure is consequently defined at the circumcenter of each tri¬ 
angular/tetrahedral element. The accuracy of the covolume scheme was estimated 
by Nicolaides 0 to be second order for a mesh with all its triangular elements 
having the same circumradii (i.e. structured-triangular mesh) and first order ac¬ 
curate otherwise. These accuracy estimates were in agreement with numerical 
experiments 0. Several forms of the covolume method were then developed for 
both two-dimensional (2D) (only on planar domains) and three-dimensional (3D) 
domains, where the difference was mainly in the convective term discretization 

mmmm. 

The conservation properties of the covolume method were later revealed by 
Perot ifTOll . The divergence form of Navier-Stokes equations was proved to con¬ 
serve the momentum and kinetic energy both locally and globally. On the other 
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hand, the rotational form of Navier-Stokes equations was found to conserve the 
circulation and kinetic energy locally and globally for both 2D iflOll and 3D ([TT| 
discretizations. These conservation properties of the covolume method, in addi¬ 
tion to the attractive properties of its differential operators that mimic the behavior 
of their continuous counterparts, shed light on the merit of using discrete calculus 
methods to solve other physics problems IIT21 . 

Another approach to develop conservative discretizations for incompressible 
flows emerged from the computer graphics community, aiming to mitigate the ef¬ 
fects of numerical viscosity that causes detrimental visual consequences [|T3l 1X4 ]. 
In this approach, the Navier-stokes equations were discretized through the dis¬ 
crete exterior calculus (DEC) framework; the discretization of the smooth exterior 
calculus operators [15, HH. A main advantage of DEC discretizations is the ap¬ 
plicability to simulate flows over curved surfaces, unlike the covolume approach. 
The resulting discrete equations had similarities with the covolume method, with 
the differences mainly in the convective term discretization. In practice, the con¬ 
vective term was not discretized using DEC but employed a method of charac¬ 
teristics with an interpolation scheme based on Kelvin’s circulation theorem lfl3l . 
or a finite volume based approach lH4ll . However, the presented numerical test 
cases using the DEC approach lacked comprehensive quantitative analysis of the 
scheme’s accuracy and its conservative behavior. 

This paper presents a discretization of the Navier-Stokes equations through 
discrete exterior calculus. Hence, similar to previous DEC-based discretizations, 
the developed discretization is capable of simulating flows over curved surfaces, 
which distinguishes it from the covolume method. In addition, the convective 
term in the presented discretization is different from previous DEC-based and co¬ 
volume discretizations. The Navier-Stokes equations are first rephrased using the 
exterior calculus notation in Section [2j The DEC discretization of NS equations 
is then derived in Section [3] for both 2D and 3D cases, highlighting its distinction 
from the covolume method and previous DEC-based discretizations. In Section 
|4j several numerical experiments are illustrated for incompressible flows over 2D 
flat/curved domains to benchmark the convergence and conservative behavior of 
the developed scheme. The paper closes with conclusions summarizing the main 
features of the presented discretization, and addressing potential future develop¬ 
ments. 
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2. Navier-Stokes equations in exterior calculus notation 

The first step in deriving a DEC discretization of NS equations is to express 
the equations using the exterior calculus notation. This is done first by starting 
from the well-known vector calculus formulation of NS equations (in Euclidean 
space) and substituting with identities relating the differential operators; viz. div, 
grad and curl, with exterior calculus operators; viz. exterior derivative, Hodge 
star and wedge product. An alternative derivation of the resultant formulation is 
then presented, starting from the coordinate invariant formulation of NS equations 
expressed in terms of the Lie and exterior derivatives. Readers not familiar with 
exterior calculus may refer to [17. 18] for a concise introduction to the topic. 

Considering the incompressible flow of a homogeneous fluid with unit density 
and no body forces, the governing equations for the fluid motion are the Navier- 
Stokes equations expressed as 


—— pAu + (u.V)u + Vp = 0, 
at 

V.u = 0, 


(la) 

(lb) 


where u is the velocity vector, p is the pressure and p is the dynamic viscosity. 
Substituting with the tensor identities: Au = V(V.u) - Vx(Vxu) and (u.V)u = 
|V(u.u) - u x (V x u), and considering the incompressibility condition V.u = 0, 
Eq. ( [Ta| ) can be expressed in its rotational form as 

dw j 

— + pV x V x u - u x ('V x u) + V// = 0, (2) 

ot 

where p d is the dynamic pressure defined as p d - p + \ (u.u). 

The notation transformation is carried out by applying the flat operator (b) to 
Eqs. ([2]) and ( |Tb| ), and substituting with the following identities 

(V x V x u) b = (-l) w+1 *d*du b , 

(u x (V x u)) b = (-l) N+l * (u b A *du b ), 

(V.u) b = *d * u b , 

(V/) b = dp d , 

where * is the Hodge star, d is the exterior derivative, and A is the wedge product 
operators. The action of the flat operator (b) on a vector u transforms it into a 
1-form u b . The sign (-1 ) v+l in the first two relationships implies a negative sign 
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only in the two-dimensional (2D) case, where N is the space dimension (i.e. N = 
2 and 3 for the 2D and 3D cases, respectively). The above relationships can be 
easily verified using the definition of differential forms and the action of exterior 
derivative, Hodge star and wedge product operators on them. Substituting with 
Eqs. ([3]) in Eqs. ([2]) and ( fib] ), the Navier-Stokes equations are then expressed as 


f) u b 

-+ (- l) N+l p * d * du b + (-l)^ 2 * (u b A *du b ) + d p d = 0, (4a) 

ot 

* d * u b = 0, (4b) 

where the velocity field is represented by the 1-form u b and p d is now the dynamic 
pressure 0-form. 

The above formulation in Eqs. Q is derived starting from the standard vector 
calculus formulation of NS equations. We now present an alternative derivation 
purely in terms of differential forms, starting from NS equations formulated using 
the exterior derivative (d) and the Lie derivative (£) operators. The Navier-Stokes 
equation in coordinate invariant form is (see [IT91 pg. 589 for Euler equation in 
this form) 


<9u b 1 

——h p(dd + dd)u b + £ u iT - -d(u b (u)) + dp = 0, (5) 

ot 2 

where 6 is the codifferential operator defined as 6 = (-1 y v,1 ” | J +l * d*, which acts 
on a k-form and results in a (k — l)-form, and N again is the space dimension. 
Therefore, the incompressibility condition (Eq. ( |4b| )) translates to du b = 0. The 
operator (dd + dd) is the Laplace operator in the exterior calculus notation, which 
differs by a negative sign from that defined in the vector calculus notation (e.g. 
the Laplace operator in Eq. (fTa|)). The Lie derivative term (£ u u b ) evaluates the 
change of the velocity 1-form u b along the velocity vector field u, and the term 
u b (u) represents the dot product of the vector field u with itself. Using Cartan 
homotopy formula (see llT9ll . pg. 430), the Lie derivative term is expressed as 

£ u u b = d/„u b + 7 u du b = d(u b (u)) + z u du b , (6) 

where i x a is the interior product of a k- form a with a vector field x. Accordingly, 
considering the incompressibility condition du b = 0, Eq. ([5]) can be expressed as 

(9ll^ 1 

—— + //ddu b + z u du b + -d(u b (u)) + dp = 0. (7) 

ot 2 
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Defining the dynamic pressure O-form as p d = p + \ (u b (u)), and substituting 
with 8 = (-1 f +1 * d* since the codifferential operator 8 in Eq. ([7]) acts on the 
2-form du b (hence k = 2), Eq. ([7]) then takes the form 


—— + (-if + f * d * du b + / U du b + d p d = 0. (8) 

at 


The interior product term i u du b can be written in terms of the Hodge star and 
wedge product using the formula (see Ifl5l Eq. 8.2.1) 


i x a = (-\f N ~ k > * (* tt a x b ). 


(9) 


for a k-form a and a vector field x. Therefore, with x = u and a = du b (hence 
k = 2), Eq. ([8]) becomes 


du* 

dt 


+ (-1 ) N+1 p * d * du° + (- 


-l) w “ 2 * (u b A *du b ) + d p d = 0, 


( 10 ) 


where the order of the wedge product was flipped using the relation a A [3 = 
(-lf/3 A a, with a = *du b is an (N - 2)-form and p = u b is a 1-form. Noting that 
(_ 1 )V+ 2 _ if -2 = (— 1)^ for any N, Eq. ( pT)] ) is exactly the same as Eq. ( |4a| ), 
which was derived earlier starting from the standard vector calculus formulation 
of Navier-Stokes equation. 

Applying the exterior derivative operator to Eq. ( [TO] ) (equivalent to taking the 
curl of the momentum equation Eq. Q), and considering the exterior derivative 
property dd = 0, the resulting governing equation is 


-+ (-if +1 ud * d * du b + (-if d * (u b A *du b ) = 0. (11) 

dt 

The DEC discretization of Navier-Stokes equations is carried out through the dis¬ 
cretization of Eq. dt The advantage of discretizing the vorticity form of the NS 
equations is highlighted in the next section. 


3. The discretization method 

In this section, the notation of the simplicial mesh used to discretize the sim¬ 
ulation domain is first introduced. This is followed by brief description of the 
discrete differential forms and some of the discrete operators. The discretization 
of NS equations is then derived for both 2D and 3D cases. The simplicial mesh 
and discrete operators concepts are discussed only briefly here. Readers interested 
in more details may refer to lH311 i~6 . 171. 
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Figure 1: A sample simplicial mesh in 2D showing the primal simplices (in black 
color) and their dual cells (in red color). The positive orientation of the primal 
2-simplices and dual 2-cells is counterclockwise. 


3.1. The domain discretization 

Let Q be the physical domain of dimension N = 2 or 3. The domain Q is 
approximated by the simplicial complex K. Following the notation in lfT51[T7Tl . a 
domain simplex cr of dimension k is denoted by cr k e K. A ^-simplex cr k is de¬ 
fined by the nodes forming it as cr k = [vo,..., v’/J. where the subscripts represent the 
nodes indices. The order of the nodes defining a simplex implies its orientation. 
The top dimensional simplices cr N are assumed to have been oriented consistently, 
whereas the orientation of the lower dimensional simplices is arbitrary. An exam¬ 
ple 2D mesh is shown in Fig. [T] The number of k -simplices in the discrete mesh 
is denoted by N k . Therefore for the mesh in Fig. [Tj A 2 = 8, N\ = 15, and A 0 = 8. 

Associated with the primal simplicial complex K is a dual complex *K con¬ 
sisting of cells. For a primal k-simplex cr k e K, its dual is an (N - k)-cell denoted 
by -kcr k G *K. The dual mesh considered here is the circumcentric dual, shown 
in Fig. |T]in red color. For a 2D mesh, the dual of a triangle is its circumcenter, 
the dual of a primal edge is the dual edge connecting the circumcenters of the 
two triangles sharing the primal edge, and the dual of a primal node is the 2-cell 
(polygon) formed by the duals of the edges connected to this primal node. For 
the case of a triangular mesh representing a curved surface, the dual edges can be 
kinked lines and the dual cells can be non planar. In 2D, the positive orientation of 
both primal 2-simplices (triangles) and dual 2-cells (polygons) is assumed to be 
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counterclockwise. The orientation of the primal 1-simplices (edges) is arbitrary, 
however their orientations induce the dual edges orientations. The dual edges can 
be oriented simply by rotating the primal edge orientation 90 degrees along the 
triangles orientation (i.e. counterclockwise), as shown in Fig. [T] In the 3D case, 
the orientations of the primal edges and faces (i.e. triangles) are arbitrary, how¬ 
ever, their orientations induce the orientations of their duals through the right hand 
rule. 

The simplicial meshes considered here are Delaunay meshes, with an extra 
requirement only for the /V-simpliccs with a face on (Kl. the domain boundary. 
Previous investigation Il20l showed that in order to correctly represent the discrete 
Hodge star operator, the mesh interior iV-simplices have to be pairwise Delaunay, 
while the iV-simplicial elements with a face on the domain boundary have to be 
one-sided (i.e., with respect to the face of the iV-simplex on the domain boundary, 
both the simplex circumcenter and its apex have to be on the same side). Such 
Delaunay meshes can easily be generated using commercial or open source mesh 
generators. 

3.2. Discrete exterior calculus 

Discrete exterior calculus provides discrete definitions to many of the exte¬ 
rior calculus operators (e.g., exterior derivative, Hodge star, wedge product, etc.) 
0211161. These discrete operators have the advantage that they satisfy many of the 
rules/identities that characterizes their smooth counterparts. Such mimetic behav¬ 
ior of the discrete operators is known to result in preserving the physics implied 
in the smooth governing equations at the discrete level Iil2l - which consequently 
improves the physical fidelity of the numerical discretization method. 

The key entities in exterior calculus are differential forms, which according to 
Flanders liTHl are best thought of as “the things which occur under integral signs” 
010. For example, considering the integration of a scalar function over a three 
dimensional space; i.e. f a dV, an example of a 3-form is Z? 3 = a dV, where 
the superscript 3 indicates a 3-form. Similar examples for 2-forms/l-forms can 
be deduced from integration of a vector field over areas/lines, whereas a 0-form 
represents a scalar function. While a smooth differential form can be integrated on 
a ^-dimensional domain, the evaluation of a discrete k-form can be thought of as 
the integration carried out on discrete k-dimensional mesh objects; i.e., line, area 
or volume. Therefore, a discrete differential k-form ultimately associates a scalar 
with a discrete k-dimensional mesh object through integration. For example, for 
the smooth velocity 1-form u', its discretization can be defined on primal edges 
cr 1 as u.dl , or on dual edges ★cr 1 as u.r/l], representing a primal or dual 
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discrete 1-form, respectively. Similarly, discrete O-forms are defined as scalars 
on the primal or dual nodes, and the discrete 2-forms are scalars resulting from 
the integration of smooth 2-forms on primal 2-simplices (triangles) or dual 2-cells 
(polygons). 

The space of discrete A-forms defined on primal and dual mesh complexes 
is denoted by C k {K) and D k (*K), respectively. Such spaces are related via the 
discrete exterior derivative and Hodge star operators as shown in the following 
diagrams for both 2D and 3D cases: 


C°(K ) C\K) — C\K) 


*0 

*1 

' J ' ' 

<• 


D 2 (+K ) ^— D ] (+K) — D°(+K) 


C°(K ) 


do 


C\K) 


di 


C\K ) 


d2 


C\K) 


*0 
\ f 

D\*K) 


*1 

*2 

^ > 

< \ 


D 2 (*K) D l (*K) v—— D°(-kK) 



( 12 ) 


(13) 


The discrete exterior derivative operator d* maps primal A-forms to primal 
(A + l)-forms. The discrete exterior derivative operator that maps dual A-forms to 
dual (k + l)-forms is the transpose of the d^-yt-i) operator, with a negative sign 
only for the dj operator in 2D (due to the defined mesh orientation convention). 
The discrete Hodge star operator * k maps primal fc-forms to dual (N - fc)-forms. 
The inverse map of the discrete * k operator is *^ 1 , which maps dual (N - A')-forms 
to primal A-forms. 

The discrete exterior derivative operator d k is a sparse matrix that is defined as 
the transpose boundary operator for the (k+ l)-simplices. For example, for the 2D 
mesh in Fig [T] the discrete di operator, that maps the primal 1-forms defined on 
the primal edges to 2-forms defined on the triangles, is an N 2 x N\ matrix defined 
as 


[diF = 


+ 1 , 
- 1 , 
0 , 


if the edge j is a face to the triangle i, 
and their orientations are consistent, 
if the edge j is a face to the triangle i, 
and their orientations are not consistent, 
if the edge j is not a face to the triangle i. 


(14) 
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For the primal nodes laying on the domain boundary; e.g. node 3 in Fig. |Tj 
the boundary of their dual 2-cells (polygons) includes primal boundary edges. 
Accordingly, the [-do] matrix, represented by the transpose boundary operator of 
these dual 2-cells, is complemented by an additional operator accounting for the 
primal boundary edges. The discrete Hodge star operator * k , on the other hand, 
is a diagonal matrix with the i-th diagonal element being the ratio between the 
volume of the dual (N - k)-simplex ★cF and the volume of its primal ^-simplex 

07; i.e. Fjp In regards to the wedge product operator, its discrete definition is 
provided within the discretization presentation in the next two subsections. 

3.3. Two dimensional discretization 

The two-dimensional DEC discretization of Navier-Stokes equations is de¬ 
rived in this subsection and the three-dimensional discretization in the next sub¬ 
section. As pointed out by Hirani et al. m, due to the intrinsic coordinate inde¬ 
pendent nature of exterior calculus, the derivation below, for the 2D case, results 
in a numerical method that works without change for both planar domains and 
curved surfaces. The dimension of the embedding space does not matter, neither 
do the details of the embedding. This characterizes a key distinction between the 
DEC-based approach and the covolume method. In the latter, the discretization of 
NS equations was commenced by taking the dot product of the momentum equa¬ 
tion and the vectors perpendicular to the triangles’ faces. Such normal vectors 
may not be unique for a simplicial mesh approximation of a curved surface. 

The discretization of NS equations is carried out here following the exact frac¬ 
tional step method HI [22]. This consists mainly of two steps, the first is to dis¬ 
cretize the vorticity formulation of NS equations (Eq. ©>. and the second is 
to substitute the velocity by its definition as the curl of a stream function, where 
the latter in its discrete manifestation becomes the unknown degrees of freedom 
in the resulting linear system. Sufficient details are provided since the numerical 
experiments presented in this paper are limited only to two-dimensional flow test 
cases. 

Discrete exterior calculus discretization of Eq. ( |TT| ) first requires the location 
on the mesh where one defines the discrete variables such that the smooth forms 
are satisfied in an integrated sense. Then, the smooth forms are replaced with their 
discrete counterparts, and the smooth operators are substituted by the appropriate 
discrete operators. Starting with the time derivative term in Eq. CD. we choose 
to define the 2-form du : ' on the dual 2-cells. Accordingly, all the other terms of 
Eq. ( [IT) are constrained to be defined on the dual 2-cells for consistency. Back 
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to the time derivative term, it then follows, according to the diagram in Eq. 0, 
that the velocity 1-form u b in this term is defined on the dual 1-cells (i.e. dual 
edges). We denote the velocity 1-form u b defined on the dual edges by u, which 
represents the integration of the velocity vector field along the dual edges; i.e. 
u = ii.c/1 e D l (irK). The velocity 1-form u may be referred to as the normal 
velocity form or the flux, since it represents the velocity normal to the triangles’ 
faces (i.e. primal edges). Similarly for the viscous term in Eq. GD. it follows 
from the diagram in Eq. ( fl2] ) that the velocity 1-form in this term is also defined 
on the dual edges. 

In regards to the convective term, because the term d * (u b A *du b ) is defined 
on the dual 2-cells, the term *(u b A *du b ) has to be defined on the dual edges, and 
therefore (u 1 ’ A *du b ) is defined on the primal edges. Since u 1 ' is a 1-form, then 
*du b is a 0-form (in 2D), and therefore the wedge product is carried out between a 
1-form and a 0-form, such that the outcome of this wedge product must be defined 
on the primal edges. This implies that both wedge product portions are defined on 
primal mesh objects. Accordingly, the velocity 1-form in the first portion of the 
wedge product term is defined on the primal edges, whereas the velocity 1-form in 
the second portion of the wedge product term is defined on the dual edges, making 
the 0-form *du b to be defined on the primal nodes. We denote the velocity 1-form 
u b defined on the primal edges by v; i.e. v = J t u.r/1 e C l (K). The velocity 1-form 
v represents the velocity tangential to the triangles edges. 

After substituting with the appropriate discrete operators, with N = 2 for the 
2D case, the discretization of Eq. ( [TTj ) takes the form 



[[-dj]£/ + d*v] = 0, (15) 


+ [-dj] W v * 


-l 

o 


where U is the vector containing the dual (normal) velocity 1-forms u for all mesh 
dual edges, and the discrete wedge product of the tangential velocity 1-form v with 
the 0-form *du is represented by the W v matrix, where the matrix W v contains the 
values of the tangential velocity 1-form v. The superscripts n and n + 1 in the time 
derivative term are due to time discretization and these subscripts are suppressed 
for the convective and viscous terms, deferring the time discretization of the vis¬ 
cous and convective terms until later in this subsection. The discrete operation 
[—do C/] evaluates the circulation of the velocity forms u along the dual 2-cells 
boundaries. Since a portion of these dual 2-cells boundaries may consists of pri- 
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mal edges, as discussed in Section [3T2] , \d h V\ then complements the velocity circu¬ 
lation, accounting for the part that depends on the velocity 1-forms v on the primal 
boundary edges. The vector V contains the tangential velocity 1-forms v for all 
mesh primal edges, and the matrix d b is then defined as d /; = 0.5|d^|c//ay(d[ 1 2 ), 
where |dg| is the matrix dj with all entries made non-negative and 1 2 is a vector 
of ones with N 2 entries, and diag{.) to be a diagonal matrix composed of the en¬ 
closed vector entries. The subscript b emphasizes the fact that d /; “closes” the dual 
2-cells for the boundary vertices by traversing along the primal boundary edges 
in the orientation direction of the dual 2-cells. Such a boundary contribution van¬ 
ishes in the time derivative term due to the time discretization. It also vanishes for 
the other dj operators in both the viscous and the convective terms, as shown in 
Appendix A 

A definition for the discrete primal-primal wedge product was developed in 
llT5Tl (pg. 74 Eq. (7.2.1)), according to which the wedge product between a discrete 
primal 1-form a and a discrete primal 0-form /3 defined over a primal edge [z, j] is 


(a A p, [; i , ;']) = - ((a, [z, j])(j3, [./']) - < a , [j, /]></?, [z']» 




(16) 


where [z, j] is the edge formed by nodes [z] and [/'], (a, [z, /']) is the discrete form 
a defined on the simplex [z, j], with the property (a, [z, j ]) = -(a, [j, z]). Recalling 
that 1 [-dj] U is a vector with N 0 entries, W v is a sparse A) x N () matrix defined as 
W v = 0.5 diag(V) |d 0 |. Accordingly, the action of this wedge product operation (in 
2D) at each primal edge is to take the average of the vorticity 1 [-dj] U evaluated 
at the edge’s nodes and multiply it by the edge’s tangential velocity 1-form v. It is 
worth noting that the vector U in Eq. ( [T?] ) includes the normal velocity 1-forms u 
for all mesh edges, including those that might be given as boundary conditions. 

In order to obtain a linear representation for the convective term, we con¬ 
sider the tangential velocity 1-forms v to be given through an interpolation of 
previously-known normal velocity 1-forms u. A velocity vector field can be cal¬ 
culated inside each triangle through the interpolation of the velocity 1-forms *^'zz 
defined on the triangle’s faces. The interpolation is carried out here using Whitney 
maps urn Since the velocity 1-forms *^u are closed forms; i.e. dj ^j 1 u = 0, 
the interpolation will result in a constant velocity vector field over each triangle 
(see Il23ll . theorem 5.4). The constant vector field obtained by Whitney map in¬ 
terpolation is the one that would result from requiring a constant vector field with 
the given fluxes of an incompressible vector field. It is worth noting that other in- 
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terpolation methods, like these in 0], can also be used to interpolate the velocity 
vector on triangles, which would result in velocity vectors not significantly differ¬ 
ent from these interpolated using Whitney maps. The tangential velocity 1-forms 
v can be then calculated on each primal edge by averaging the velocity vector 
fields on the triangles sharing the edge. 

The second step in the exact fractional step method is to substitute for the ve¬ 
locity 1-forms u (in Eq. ( fl5| )) by its definition as the curl of a stream function. The 
incompressibility condition (Eq. ( |4b| )) in DEC notation is expressed as * 2 d| *(' U. 
This implies that the vector U is in the null space of the matrix [* 2 d| *)"' ]■ Since 
[* 2 di*7 1 ][*id 0 ] = * 2 did 0 = 0, the columns of the matrix [*id 0 ] then contain a 
basis of the null space of [* 2 dj*i *]. Therefore, the vector U can uniquely be ex¬ 
pressed in terms of the basis [*id 0 ]; i.e., U = *id 0 T / . In vector calculus notation, 
this is equivalent to expressing a divergence-free velocity vector as the curl of a 
stream function. According to the diagram in Eq. ( [T2| ), the vector V F includes the 
stream function if/ defined as 0-forms on the primal mesh nodes (i.e. if/ e C°(K)). 
Substituting with this representation of U, Eq. ( fT5] ) becomes 


^[-d 0 r ] *! d 0 T" ,+1 -n[- dj] *, d 0 *o 1 [-do] *i do'F 

+ [-d 0 r ] *i W v *o' [-d T 0 ] do'F = F, (17) 

with the vector F = ^[-dj]t/" + /u[- dj] *i d 0 *,]' d b V - [-dj] *i W v d b V. Eq. 
( [T7| ) describes the resultant linear system to be solved. The degrees of freedom 
in the above linear system are the stream function 0-forms (scalars) defined on 
the primal mesh nodes. Therefore, the resulting system is a sparse Nq x Nq linear 
system. 

For the current 2D case, it is worth noting the correspondence between the 
velocity 1-form u and the mass flux across the primal edges. While the discrete 
velocity form u is formally defined as the integration of the velocity field on the 
dual edges (u = f i u.dl), it can be used to approximate the mass flux normal 
to the primal faces (edges) as Uf = u. The incompressibility condition then 
implies a zero summation of the mass fluxes across the faces (edges) of each 
triangular element; i.e., d\Uf = dif*]" 1 ^] = 0. On the other hand, if the mass flux 
across the primal edges is the known quantity, it can be used to approximate the 
velocity 1-form u as u - *\Uf. Recalling that the actual degrees of freedom in 
the present discretization are the stream functions if/ and the velocity 1-form u is 
calculated as u = *|d 0 (A, the mass flux across the primal edges is then uj = do if/. 
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This implies that it is the mass flux Uf that is calculated first and is then used to 
approximate the velocity 1-form u as u = *\Uf = *id 0 t/c Such a correspondence 
between the velocity 1-form u and the mass flux uj provides some flexibility in the 
case that an initial analytical expression for the velocity vector field defined on a 
smooth surface is available. In order then to set the initial solution on the discrete 
triangulation objects, it is possible to calculate Uf by integrating the mass flux 
normal to hypothetical curved primal edges and then approximate u as u = *\Uf. 
Otherwise, the velocity 1-form u can be initialized by integrating the velocity field 
on hypothetical curved dual edges. It is worth noting that the integration of the 
mass flux might be preferred in order to accurately guarantee zero net mass flux 
across the domain boundaries. 

We now elaborate on the time discretization wherein the linear system in Eq. 
( fT7| ) is solved in two steps as a predictor-corrector method. First, we advance the 
system explicitly by a half time step 


0.5At 


[-dn] *1 do 


vp»+l = p + 


H[-dl] *! d 0 V [—dj] 


-[-dj] *1 Vh; 


v *0 


1 [-do r ] 


<id 0 T"\ (18) 


where the matrix IT" incorporates the tangential velocity forms v at time n. Af¬ 
ter solving the linear system ( fT8| ), the normal velocity 1-forms at time n + \ are 
calculated as U n+1 = *id 0 v F' !+ T These normal velocity 1-forms are then used to 
predict, at time step n + \, the velocity vector field at each triangle element through 
Whitney map interpolation [fTTl . The tangential velocity 1-form v is then calcu¬ 
lated at each primal edge as the average of the velocity vectors in the neighboring 
triangles, in the direction of the primal edge, multiplied by the primal edge length. 

This results in the tangential velocity 1-forms v at time n + and then W v 2 . The 
latter matrix is then used to calculate the new velocity 1-forms u at time n + 1. The 
evaluation of the tangential velocity 1-forms v at half time step (n + is expected 
to improve the kinetic energy conservation, as was shown before by Perot [HOI for 
the covolume method. The second linear system to solve is then 


At 


[-dj] *i d 0 -yu[-do] *, d 0 * 0 * [-dj] *j d 0 


+[-dj] WT 5 *o 1 [-dj] *i d Q 


T" ,+1 = F, (19) 
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where the right hand side vector F in both Eqs. ( [T8] ) and ( p~9| ) also contains the 
contribution from any stream function 'F boundary conditions. 

The solution of the Navier-Stokes equations through the exact fractional step 
method is known to significantly reduce the size of the solved linear system 
0[|22|. While the discretization of the momentum equation ( fTO] ) would result 
in a number of degrees of freedom equals to the number of edges plus the number 
of triangles; i.e. (N\ + AT), the presented discretization has a number of degrees of 
freedom equals only to the number of primal nodes N (] . In addition to reducing the 
linear system size, the presented discretization always maintains the incompress¬ 
ibility error at the machine precision. Recalling the discrete incompressibility 
condition * 2 cli U = * 2 di ♦ido'P = * 2 did 0 v P, and since dido = 0, the result¬ 
ing formulation guarantees the mass conservation up to the machine precision, 
regardless of the error incurred during the linear system solution. On the other 
hand, the solution of the Navier-Stokes equations in terms of a stream function 
might simplify or complicate the implementation of some boundary conditions. 
A discussion regarding the implementation of various boundary conditions, in¬ 
cluding the boundary conditions for interior walls (domain boundaries with zero 
in/out flow), through the stream function can be found in |f22| . For the case of 
interior boundaries with nonzero in/out flow, the boundary conditions can be im¬ 
plemented through Hodge decomposition manipulation for the stream functions. 
More details regarding the implementation of such boundary conditions might be 
addressed in future publications. 

The presented discretization of NS equations has similarities with previous 
discretizations. The viscous term discretization is similar to that previously de¬ 
veloped through the covolume method 0, :b. 8. j9, H)| and DEC-based [fT3l [14] 
discretizations. However, the discretization of the convective term via the inte¬ 
rior product definition in Eq. <(9]) and the wedge product definition in Eq. < [T6| ) 
makes the present discretization different from previous DEC-based discretiza¬ 
tions that adopted Lie derivative advection techniques ||T3|| or a finite-volume- 
based approach Ifl4ll in discretizing the advection term. It also distinguishes the 
present discretization from most of the covolume method discretizations, with 
some similarity with the covolume discretization developed by Perot [fill only in 
the special case of structured-triangular meshes; i.e. when the center points of the 
primal edges and their dual edges coincide. Nevertheless, unlike all covolume dis¬ 
cretizations, the current discretization is capable of simulating flows over both flat 
and curved surfaces. In addition, the current manipulation of the convective term, 
through the discrete wedge product operator, gives insight into other important 
research themes. For example, it paves the way to exploring the discretization 
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of Navier-Stokes equations through the finite element exterior calculus method 
[23:, [24]. Furthermore, it gives insight into the discretization of other convective 
terms in different physics problems; e.g. the magnetohydrodynamics governing 
equations. 


3.4. Three dimensional discretization 

The three-dimensional DEC discretization of NS equations is briefly presented 
in this subsection. In the 3D situation, the primal mesh consists of 3-simplices 
(tetrahedra), 2-simplices (triangles), 1-simplices (edges) and O-simplices (nodes). 
The duals to these primal mesh entities consist of 0-cells (dual nodes), 1-cells 
(dual edges), 2-cells (polygons) and 3-cells (polyhedra), respectively. The space 
of the discrete k-forms defined on primal/dual mesh complexes is defined accord¬ 
ing to the diagram in Eq. ( [T3| ). 


Following the same methodology in Section 3.3, we start the discretization 
of Eq. ( fTT) by choosing to define the 2-forms du’, in the time derivative term, 
on the dual 2-cells (polygons). It then follows that the velocity 1-forms u' 1 in 
this term are defined on the dual 1-cells (i.e. dual edges). The velocity 1-forms 
defined on the dual edges are denoted by a. It also follows, based on the diagram 
in Eq. ( [T3] ), that the velocity 1-forms in the viscous term are also defined on the 
dual edges. In regards to the convective term, since d * (u b A *du b ) is defined 
on the dual 2-cells, then *(u b A *du b ) has to be defined on the dual 1-cells (dual 
edges), and therefore (iT A *du b ) is defined on the primal 2-simplices (triangles). 
Since u“ is a 1-form, then *du b is also 1-form (in the 3D case), and therefore 
the wedge product is carried out between two 1-forms, where the outcome of this 
wedge product is defined on the primal triangles. This implies that both the wedge 
product portions are defined on the primal edges. Accordingly, the 1-form *du b is 
defined on the primal edges, which makes u" in this portion to be defined on the 
dual edges. Finally, the first velocity 1-form in the wedge product term is defined 
on the primal edges. We denote the velocity 1-forms defined on the primal edges 
by v. 

Substituting with the appropriate discrete operators, with N = 3 for the 3D 
case, the discretization of Eq. (jTT|) then takes the form 


d[ U Af U + yud[ * 2 di d\U - d[ *2 W v *~ l d[f/ = 0, (20) 

where U is a vector containing the dual (normal) velocity 1-forms u for all mesh 
dual edges, and W v is a sparse matrix representing the action of the discrete wedge 
product operator. It is worth noting that similar to Eq ( p~5T ), the operator dj needs 
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to be complemented by another operator that closes some of the dual 2-cells on 
the domain boundary. Such a complementary matrix is omitted in Eq. ( |20| ) for 
simplicity, but its existence should always be considered during numerical imple¬ 
mentations. 

According to the definition in [IT51 (pg. 74 Eq. (7.2.1)), the wedge product of 
a discrete primal 1-form a and a discrete primal 1-form /3 defined over a primal 
2-simplex (triangle) [i, j, k\ is 


(a A p, [i, j, k]) = a , [i, j\m [j, k]) - (a, [i, k])(J3, [k, j ]) 
b 

-(a, [j, i])(P, [i, k]) + (a, [j, k])(j3, [k, /]> 
+<a, [k, i])(j3, [i, j ]> - (a, [k, ;']></?, [j, /])), 


( 21 ) 


with the property (a, [i, j]) = -(a, [j, /]). Recalling that is a vector of 

N\ rows, W v is then a sparse N 2 x N\ matrix. For each row in W v correspond¬ 
ing to a primal triangle, the only non-zero entries in this row are for the primal 
edges belonging to this triangle. For a primal triangle [i,j,k], the matrix entry 
corresponding to the edge [/,./] is |«v, [ k , i ]) - <v, [j, k ])). 


Following the procedure in section 3.3, the velocity 1-forms are substituted 


by their definition in terms of a stream function; i.e. U = * 2 di v P. According to 
the diagram in Eq. ( [T3] ), the vector 'P includes the stream functions if/ defined as 
1-forms on the primal mesh edges (if/ 6 C 1 (70). Substituting in Eq. ( [20] ), the 
resulting linear system is 


d[^ *2d il P' ,+1 +/id[ * 2 di*i 1 d[ * 2 d i v F-df * 2 W r v * 1 1 d[ = d[ — U n . (22) 


The degrees of freedom in the above linear system are the stream function 
1-forms (scalars) defined on the primal edges. Therefore, the resulting system is 
a sparse N\ x N\ linear system. The time discretization of Eq. ( |22[ ) can then be 
implemented similar to the 2D case. 

The presented 3D discretization of the viscous term in Eq. ( |22[ ) is similar 
to previous DEC-based Ifl3l 14] and covolume method []7l 111] discretizations. 
However, using the interior product definition in Eq. ([9]) and the discrete wedge 
product definition in Eq. ( |2Tj ) makes the present discretization different from all 
previously developed 3D DEC-based/covolume discretizations. 
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4. Results and Discussion 


In order to benchmark the performance of the presented discretization, sev¬ 
eral 2D simulation experiments are performed for flows over both flat and curved 
surfaces. During all simulations, Eqs. ( |T8| ) and ([T9]) are solved consecutively at 
each time step, where direct LU decomposition solver is used to solve the linear 
systems. As pointed out earlier, the mass conservation is guaranteed by the dis¬ 
cretization construction. Vorticity is also locally conserved due to the discretiza¬ 
tion construction, as shown earlier by Perot IfTOll . Global vorticity conservation up 
to the machine precision was observed during all conducted simulations, where in 
the presence of solid walls the vorticity flux comes only from the no slip bound¬ 
aries as it should be for incompressible flows. Therefore, the results presented 
below mainly quantify the numerical convergence rate of the discretization and 
the conservation of the kinetic energy. 


4.1. Driven cavity 

Driven cavity simulations are carried out on a unit square domain at Reynolds 
number (Re) of 1000. Solid wall boundary conditions are imposed on the left, 
right and bottom boundaries. The top boundary has zero flux (i.e. u), and a unit 
tangential velocity (i.e. v) boundary conditions. Therefore, the stream functions 
on all boundary nodes are set to an arbitrary constant value. The fluid dynamic 
viscosity (ju = \/Re in our normalized units) is set to 0.001, and the time step is 
At = 0.1. The simulations are carried out on a Delaunay mesh and a structured- 
triangular mesh (consisting of isosceles right triangles) with 32482 and 32258 ele¬ 
ments, respectively, which has almost the same resolution as a 128 x 128 Cartesian 
mesh. 

Fig. [2] shows cross-sections at the domain center lines for the steady state 
velocity profile at simulation time T = 100. The results are compared with well- 
established simulations by Ghia et al. fl25ll for Re = 1000 using a 128 x 128 
Cartesian mesh. The comparison shows an agreement with Ghia’s results for both 
mesh types, which reflects the numerical solution fidelity. 


4.2. Taylor-Green vortices 

The simulation of Taylor-Green vortices is carried out on a square domain of 
dimension \—n, n\ in both x and y directions. The decay of Taylor-Green vortices 
with time has an analytical solution that for the 2D case is expressed as lf26l [271 

u x = - cos(x) sin (y)e~ 2vt 
u y = sin(x) cos (y)e~ 2vt 
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Figure 2: Cross-section of the velocity profile at the two domain center lines for 
driven cavity test case at Reynolds number = 1000. 


with v to be the kinematic viscosity. The simulation is conducted using a Delaunay 
mesh consisting of 50852 elements with periodic boundary conditions applied on 
all domain boundaries. This requires only to fix the stream function at one primal 
node to an arbitrary value in order to get a unique solution. The simulation is 
carried out using a time step At = 0.1 and kinematic viscosity v = 0.01. 

Fig. [3a]shows the vorticity contour plot for Taylor-Green vortices at simulation 
time T = 10. A cross section of the velocity y-component u y along the horizontal 
domain center line is shown in Fig. 3b The simulation velocity profile is in good 
agreement with the analytical solution, as shown in Fig. [3bJ This represents a 


qualitative indication of the reliability of the current numerical implementation to 
reproduce the evolution of such unsteady flow with time. 


4.3. Poiseuille flow 

Poiseuille flow simulations are carried out to investigate the numerical con¬ 
vergence rate of the developed discretization. The simulations are conducted on 
a unit square domain. Solid wall boundary conditions are imposed on the top 
and bottom boundaries, while parabolic in/out flow conditions are imposed on the 
left/right boundaries. Therefore, by fixing the stream function at one boundary 
node, the stream functions at the rest of the boundary nodes can be directly cal¬ 
culated based on the in/out flux (i.e. u) boundary conditions. The simulation is 
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(a) 



Figure 3: ([a]) The vorticity contour plot for Taylor-Green vortices at time T = 10. 
(|b]) Cross-section of the velocity y-component profile at the horizontal center line 
for Taylor-Green vortices at time T = 10. 


carried out for structured-triangular, Delaunay and well-centered meshes of dif¬ 
ferent resolutions. The well-centered mesh is a Delaunay mesh that is optimized 
to make the circumcenter of each triangle to reside inside the triangle itself ll28l . 
The exact solution of the velocity vector field is given by u = [y(l - y),0]. 

The L 2 -norm of the velocity 1-form (u) error (see Hall el al. ®|) is calculated 

2 . F/2 


as \\u 


- u\\ 


/ u exact_ u ^ 
\ |crl| ) 


CT 


★ CT 


, and its convergence with the mesh 


elements size is shown in Fig. |4j It is observed that the velocity 1-form (flux) 
error converges with a second order rate for the structured-triangular mesh case, 
and with a first order rate otherwise. This is in agreement with previous theoret¬ 
ical analysis by Nicolaides []5l for the covolume method. Such analyses showed 
that a necessary condition to obtain a second order convergence rate is to have 
the midpoint of each primal edge to coincide with the midpoint of its dual edge, 
which is satisfied only for a structured-triangular mesh or a mesh consisting of 
equilateral triangles. The observed convergence rates are therefore in agreement 
with theory. Regarding the unstructured meshes, the well-centered mesh error is 
slightly smaller than the Delaunay mesh error, where both the well-centered and 
Delaunay mesh implementations converge in a first order fashion. 

The convergence of the interpolated velocity vector field is also investigated. 
The velocity vector field is calculated inside each triangle through the interpola- 
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Figure 4: The numerical convergence of the velocity 1-form (flux) and the inter¬ 
polated velocity vector for the Poiseuille flow test case. The dashed lines represent 
the 1-st and 2-nd order slopes. 
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tion of the velocity 1-forms l u defined on the triangle’s faces using Whitney 
maps E2. As pointed out earlier in section |3.3| such interpolation results in a 
constant velocity vector field over each triangle, implying a first order interpola¬ 
tion scheme. The velocity field is interpolated over all triangles and the L 2 -norm 
of the velocity vector error is calculated. The convergence of the velocity vector 
error with the mesh size is shown in Fig. |4j A first order convergence rate is 
observed for all considered mesh types. Although the structured-triangular mesh 
exhibited a second order convergence for the flux 1-forms, the interpolated veloc¬ 
ity vector converges with a first order rate. This can be attributed to the first order 
velocity interpolation scheme, which seems to dominate the velocity vector error. 
This is confirmed by calculating the velocity vector through the interpolation of 
the exact flux 1-forms (calculated by integrating the velocity analytical solution 
over the dual edges), which also converges with a first order rate, as shown in Fig. 

m 


4.4. Double periodic shear layer 

The simulation of a double periodic shear layer is carried out for an inviscid 
flow (p = 0) over a square domain of unit edge length. The initial flow repre- 
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sents a shear layer of finite thickness with a small magnitude of vertical velocity 
perturbation. The initial velocity vector field is expressed as Il29l 


Jtanh((y - 0.25 )/p), fory < 0.5, 

|tanh((0.75 - y)/p), fory >0.5, (24) 

u y = 8 sin(2 nx). 


with p = 1/30 and 8 = 0.05. The initial velocity 1-forms u are approximated by 
integrating mass flux normal to the primal edges and then multiplying this flux by 
the discrete Hodge star operator *|. Periodic boundary conditions are imposed on 
all domain boundaries. Therefore, it is only required to fix the stream function at 
one primal node to an arbitrary value in order to get a unique solution. 

Five simulations are conducted using a time step of At = 0.001 on structured- 
triangular meshes with number of elements equal to 3042, 12482, 32258, 50562 
and 204800. Fig. [5] shows the evolution of the vorticity contour plot with time, 
using the finest mesh. At time T = 0.8, two vortices appear to be well resolved. 
The shear layer connecting the coherent vortices becomes thinner with time and 
within a finite time interval reach the resolution of the mesh after which dispersion 
error is manifested as mesh level oscillations. The vorticity contour plot in Fig. [5] 
exhibit similarities with previous simulations by Bell et al. |[29|| . 

The convergence of the kinetic energy error with the mesh size is investigated. 
The kinetic energy is calculated as f Q u.u dQ, where the integration is carried out 
over the entire simulation domain. The velocity vector is calculated in each tri¬ 
angular element via Whitney map interpolation, as described before. The kinetic 
energy relative error ( g£f ^~ ( ^ (r) ) is then calculated at time T = 2.0 and plotted 
versus the mesh characteristic length in Fig. [5d} Except for the coarsest mesh 


case, the kinetic energy relative error converges in a second order fashion with the 
mesh size, which is expected from a scheme that is second order for structured- 
triangular meshes. Overall, the kinetic energy relative error is modest, with a 
0.3% error for the coarsest mesh (equivalent to a 40 x 40 Cartesian mesh) and 
only 0.01% error for the finest mesh (equivalent to a 320 x 320 Cartesian mesh). 
For the mesh with 50562 triangular elements (equivalent to a 128 x 128 Cartesian 
mesh), the kinetic energy relative error is 0.039%, almost one order of magnitude 
lower than a second order collocated mesh scheme using almost the same mesh 
size If29t 
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Mesh characteristic length 


(c) (d) 

Figure 5: The vorticity contour plot for double periodic shear layer with a mesh 
of 204800 elements at time: @ T=0.0, (|b| T=0.8 and @T=1.2. (|dj) The conver¬ 
gence of the relative kinetic energy error ( KE(0 f ! E ^ {T) ) with the characteristic mesh 
length at simulation time T = 2.0. The dashed red line represents the second order 
slope. 
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4.5. Taylor vortices on flat surfaces 

Two Taylor vortices are simulated for an inviscid flow (ju = 0) over a flat 
square domain of dimension \-n, n] in both directions. The vorticity distribution 
for each vortex is expressed as [|3Q|| 


oj(x,y) = — (2- y exp(o.5(l 


a \ a 2 \ \ a 2 


(25) 


with G = 1.0, a = 0.3 and r is the distance between any field point and the vortex 


center. The vorticity distribution in Eq. ( [25] ) ensures that the net circulation of 
each Taylor vortex is zero. 

The domain is initialized with a vorticity distribution for two vortices sep¬ 
arated by a distance of 0.8. Such a separation distance is just above the critical 
bifurcation distance, below which the two vortices would merge OTTi fTi I. The vor¬ 


ticity values are assigned to the primal nodes according to Eq. ( |25] ). The velocity 
1-forms u are determined by solving the Poisson equation 


* 


- ] d r 
0 u o 


: i do'P = X, 


(26) 


where X is the vector containing the known vorticity values, and T is the vector 
containing the unknown stream functions on the primal nodes. No-flux Dirichlet 
boundary conditions are imposed on the domain boundaries during the Poisson 
equation solution. 

The Poisson equation is solved only once initially and the velocity 1-forms are 
then calculated as U - *|d 0 v F. Such velocity 1-forms are used as the initial state 
for the simulations. When simulating the evolution of the two Taylor vortices, 
periodic boundary conditions are imposed on all domain boundaries. Therefore, 
it is only required to fix the stream function at one primal node to an arbitrary 
value in order to get a unique solution. 

The simulations are carried out on a mesh consisting of 132204 equilateral 
triangles, using various time steps in the range [1.0 - 0.002]. Fig. [6] shows the 
vorticity contour plot evolution with time, using a time step of 0.005. The two 
vortices initially approach and turn over each other. The vortices then move apart, 
as expected, with a thin vortex sheet connecting them that disappears at longer 
simulation time. 

The relative kinetic energy error is calculated at simulation time T = 20.0 and 


is plotted versus the time step in Fig. 6d The figure shows a second order conver¬ 
gence of the relative kinetic energy error over the entire range of time steps. Such 
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Figure 6: The vorticity contour plot for two Taylor vortices using a mesh consist¬ 
ing of 132204 elements and a time step At = 0.005 at time: 0 T=0.0. 0 T=3.0 
and 0 T=5.0. ([d]) The convergence of the relative kinetic energy error with the 
time step at simulation time T = 20.0. The dashed red line represents the second 
order slope. 0 The evolution of the relative kinetic energy with time for two 
Taylor vortices on a hexagonal domain. 
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time convergence rate can be due to the two-step time stepping scheme imple¬ 
mented during the current discretization. For practical time steps that can resolve 
the physics of the considered problem (e.g. At < 0.01), trivial relative kinetic 
energy error, below 0.01%, is observed. 

Another simulation is carried out for two Taylor vortices on a hexagonal do¬ 
main, aiming to compare with the Eulerian formulation developed by Pavlov et 
al. OH . The hexagonal domain has a side length of n and is meshed with 55296 
equilateral triangular elements. The initial condition for the two Taylor vortices 
is the same as in Eq. ( [25] ) (with G = 1.0 and a = 0.3), and periodic boundary 
conditions are applied on all domain boundaries. Since the authors in PH did 
not state the time step they used, we choose to use a time step of 0.01 during this 
test case. It is expected that the time step used in PT1 will not be significantly 
larger than 0.01 since a simulation with a time step of 0.02 was tested, however 
the time step was not small enough to capture the expected behavior of the two 
vortices, and the two vortices merged together when using the 0.02 time step. Fig. 


6e shows the change in the relative kinetic energy for extended simulation time 


up to T = 100. It is observed that the kinetic energy error is within the range of 
0.003%, which is almost three orders of magnitude smaller than the error of 2.0% 
reported in Pavlov et al. PH . 

was 


It is worth noting that the initial increase in the kinetic energy in Fig. 6e 


also observed for the previous test cases of Taylor vortices on a square domain, but 
was comparable to the subsequent decrease in the kinetic energy only for the time 
steps < 0.01. Therefore, for these time steps < 0.01, the relative kinetic energy 

was calculated by comparing with the kinetic energy 


error (™^) in Fig. 
peak at time ~ 1.0 insteacTof the initial kinetic energy KE( 0). This should better 
represent how much kinetic energy loss the implemented scheme allows. 


6d 


4.6. Taylor vortices on a spherical surface 

The ability of the current discretization to simulate flows over curved surfaces 
is explored for an inviscid flow test case over a spherical surface. The spherical 
surface, with radius equal to 1.0, is approximated by a simplicial mesh consisting 
of flat triangles connecting the primal nodes, where the primal nodes are posi¬ 
tioned on the smooth spherical surface (at radius = 1.0). For each primal edge, 
its dual is the kinked line connecting the circumcenters of the two neighboring 
triangles through the primal edge’s midpoint. During the discrete Hodge star op¬ 
erator calculations, the length of a dual edge is its length as a ki nked line. On the 
other hand, for each primal node, its dual is the non-planar polygon consisting of 
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sector contributions from all the flat triangles sharing the primal node, and its area 
is calculated accordingly as a non-planar surface area. 

The test case investigated in this section is for two Taylor vortices initially 
positioned on a spherical surface and their evolution with time is simulated. The 
spherical surface domain is initialized with two vortices each have the distribution 
given in Eq. ( [25] ), with G = 0.5, a = 0.1. When calculating the vorticity at any 
mesh node via Eq. ([25]), the distance r is measured along the sphere surface; i.e. 
geographical distance. The centers of the two vortices are separated by a distance 
of 0.4. The simulation is carried out using a mesh containing 327680 triangular 
elements. 

In order to recover the velocity 1-forms from the vorticity distribution, the 
Poisson equation ( |26[ ) is solved. During the Poisson equation solution, the stream 
function at one primal node need to be fixed in order to obtain a unique solution. 
Using the resulting velocity 1-forms as an initial condition, the evolution of the 
two vortices is then simulated using various time steps in the range [1.0 - 0.05]. 

Fig. [7] shows the evolution of the vorticity contour plot with time. Again, the 
two vortices move apart with a thin vortex sheet connecting them. The conver¬ 
gence of the kinetic energy relative error with the time step is investigated after 
simulation time T = 10.0, as shown in figure [TdJ Similar to the flow over a flat 
surface, a second order rate, on average, is observed for the convergence of the ki¬ 
netic energy relative error with the time step. This again can be due to the two-step 
time stepping scheme implemented during the current discretization. 


4.7. A ring of vortices on a spherical surface 

The behavior of a ring of N equidistant point vortices, having the same strength, 
positioned on a circle with fixed latitude on a spherical surface was previously in¬ 
vestigated theoretically 071 . It was shown that with such a configuration, the 
vortices will rotate around the z-axis in a stable fashion given that the circle’s lat¬ 
itude 6 is below a critical value and the number of vortices N < 7. For N = 6, the 
critical polar angle 6 C ~ 0.464 [f33ll . The behavior of such a ring of vortices is sim¬ 
ulated, where the point vortices are replaced with vortices having the distribution 


co = 


T 

cosh 2 (—) 

v a 7 


(27) 


with t to be the vortex strength, a is the vortex radius, and r is the distance between 
any field point and the vortex center. 

Six identical vortices, having a strength r = 3 and a radius a - 0.15, are 
placed on a unit sphere at latitude 6 = 0.4. In order to satisfy the condition 
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(C) (d) 

Figure 7: The vorticity contour plot for two Taylor vortices on a spherical surface 
meshed with 327680 elements at time: (§) T=0.0, (§> T=2.0 and (§ T=10.0. © 
The convergence of the relative kinetic energy error with the time step at simula¬ 
tion time T = 10.0. The dashed red line represents the second order slope. 
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that the integration of the vorticity over a spherical surface is zero, an additional 
vortex, with a strength r = -18 and a radius a = 0.15, is placed at the south pole 
(0 = n). The spherical surface is meshed with 81920 elements, and the simulation 
is conducted for an inviscid flow with a time step At = 0.005. 

Figs. 


8a and 8b show the vorticity contour plots at time T = 0 and T 


36, respectively. It is observed that the vortices positions seem unchanged after 
such simulation time, with some flow fluctuations around the vortices due to the 
inviscid nature of the flow. The cyclic rotation of the vortices around the z-axis 
can be detected by monitoring the relative solution change, with respect to the 
original solution, with time. Recalling that the vector U{t) contains the fluxes over 
all mesh edges at time t, the relative solution change is then defined as ll ^||^|| Q)l1 • 
Such relative solution change should be equal to zero each time the six vortices 
rotate by an angle n/3 around the z-axis. The relative solution change versus time 
is shown in Fig. [8cj which reveals the periodic nature of the vortices motion. The 
six vortices perform a n/3 rotation around the z-axis in a time period of almost 12 
time units. Accordingly, at time T = 36, the six vortices have rotated by an angle 
n around the z-axis. A small non-vanishing relative solution change, of almost 
0.01, is observed after each cycle, which is due to the developing flow fluctuations 
around the vortices, as was shown in Fig. 8b Finally, the vorticity strength along 
the circle with latitude 0 = 0.4 is shown in Fig. 8d at simulation times T - 0 and 
T = 36. The figure indicates more quantitatively that at time T - 36 the vortices 
positions are similar to the original positions due to their n rotation around the 
z-axis. The vorticity strength drop at the center of all the six vortices is due to 
the developed flow fluctuations around the vortices. Recalling that the vorticity 
integration over the spherical surface is always by construction maintained at zero, 
and noting that the strength of the single vortex at the south pole only changed by 
0.002% at time T = 36, the vorticity developed around the six vortices due to 
flow fluctuations is compensated from the six vortices themselves. In regards to 
the kinetic energy, the relative change in the kinetic energy at time T = 36 is 

KE(T=0)-KE(T=36) _ Q A Y i n-6 

ke(t= 0 ) y.u x ru . 


5. Conclusions 

A conservative discrete exterior calculus discretization of Navier-Stokes equa¬ 
tions was developed. The Navier-Stokes equations were first rewritten using the 
exterior calculus notation. The expression of Navier-Stokes equations using the 
exterior calculus notation was derived from the standard vector calculus formu¬ 
lation and verified against the coordinate invariant form in terms the exterior and 
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(b) 




Figure 8: The vorticity contour plot for 6 vortices on a spherical surface at latitude 


6 = 0.4 at time: © T =0-0 and (© T=36.0. © The relative solution change 
( "^i Q> " ) versus the simulation time, (d) The vorticity strength along a circle at 


latitude 6 = 0.4. 
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Lie derivatives. 

The discretization was carried out through the substitution by the discrete exte¬ 
rior calculus operators defined on simplicial meshes. Both 2D and 3D discretiza¬ 
tions were developed, but with numerical implementations only for flows over 
surfaces. The main distinction between the developed discretization and previ¬ 
ous unstructured conservative discretizations was in the convective term. In the 
2D case, the current convective term discretization is different from all previous 
DEC-based discretizations and most of the covolume method discretizations. The 
developed discretization has similarities with the covolume discretization by Perot 
[flOll only for the special case of a structured-triangular mesh. Nevertheless, unlike 
all covolume discretizations, the presented discretization is capable of simulating 
flows over both flat and curved surfaces. In regards to the 3D case, the current 
convective term discretization is different from all previous unstructured conser¬ 
vative discretizations. An additional merit of the presented methodology is the 
manipulation of the convective term through algebraic discretization of the in¬ 
terior product operator and a combinatorial discretization of the wedge product. 
Such approach paves the way to explore the application of the finite element ex¬ 
terior calculus method to discretize Navier-Stokes equations. Moreover, it gives 
insight into the discretization of similar convective terms in other physics prob¬ 
lems; e.g. the magnetohydrodynamics governing equations. 

Several 2D simulation experiments were carried out to benchmark the dis¬ 
cretization performance. The convergence of the velocity 1-forms (i.e. fluxes) 
was found to be of second order for structured-triangular meshes and of first order 
otherwise. This is in agreement with previous theoretical estimations developed 
for the covolume method. In regards to the conservation properties, due to the dis¬ 
cretization construction, both the mass and the vorticity are conserved locally and 
globally up to the machine precision. The kinetic energy relative error converged 
in a second order fashion with the mesh size for flows over flat surfaces. The 
convergence of the kinetic energy relative error with the time step was also found 
to be of a second order for flows over both flat and curved surfaces. Such conser¬ 
vation properties, the ability to simulate flows over both flat and curved surfaces 
and the relatively small size of the linear system make the presented discretization 
attractive for both physics and engineering applications. 
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Appendix A. The complementary contribution to the dual 2-cells boundary 
operator 


After defining the velocity 1-forms on primal/dual mesh entities and substitut¬ 
ing with the appropriate discrete operator, the discretized Navier-Stokes equation 
was expressed, as in Eq. ([15]), as 


TJ n +1 _ JJ n 

[~<]-—-M~dp] *i d 0 *o 1 [[—dj]C/ + d b v] 

+ [-<£] *i W v *~ 0 l [[-<%]£/ + d b v] = 0, (A.l) 

For the dual 2-cells touching the domain boundary; e.g. the 2-cell whose 
dual is the primal nodes 0, 1, 3, 4, .. in Fig. [I] the discrete operators [-dj] 
in the viscous and convective terms are complemented by the operator d b that 
closes such dual 2-cells boundaries by the primal boundary edges. Such domain 
boundary contribution vanishes, however, for other [-dj] operators, underlined in 
the above equation. Regarding the time derivative term, since the entries of the 
tangential velocity forms vector V are calculated at an intermediate time step, the 
domain boundary contributions complementing [-dg]t/" +1 and (-d^ ] IJ" will then 
cancel each other. 

In regards to the viscous term, starting from its smooth exterior calculus form; 
i.e. //d * d * du b , it can be expressed as /uda, with a = *d * du b . Considering 
the domain boundary contribution, the discrete viscous term is then expressed 
as /i [-dj]A + d /; A'], with A to be the vector containing the discrete a 1-forms 
defined on the dual edges, and A' as the vector containing the discrete 1-forms 
a'. Similar to a, the smooth form a' is defined as a' = *d * du b , whereas its 
discrete version is defined however on the primal edges. It follows accordingly, 
based on the diagram in Eq. ( fT2] ), that the discretization of du b included in the 
a' form is defined on the primal triangles. Since the smooth velocity form u b is 
retrieved through Whitney map interpolation as a constant form over each triangle, 
the exterior derivative of such constant velocity form then vanishes; i.e. J r du b = 
0. This implies that a', and therefore the domain boundary contribution, vanishes 
for the viscous term. Using a similar argument for the convective term d * (u b A 
*du b ), it follows that the discretization of du b , in the wedge product, is also defined 
on the primal triangles, and therefore is equal to zero. Accordingly, the domain 
boundary contribution to the first (underlined) [-dg] operator in the convective 
term also vanishes. 
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